! This case demonstrates the effects of TVD modifications on Lax-Wendroff scheme.
!
! Weighted-averaged flux:
!
!   f_{i+1/2} = (1 + 2 * a * c) / 2 * (u * rho_i) + (1 - 2 * a * c) / 2 * (u * rho_{i+1}).
!
! When a = 1/2, it is thhe Lax-Wendroff flux:
!
!   f_{i+1/2} = (1 + c) / 2 * (u * rho_i) + (1 - c) / 2 * (u * rho_{i+1})
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2022-03-13: Initial creation.

program tvd_adv_1d_case

  use adv_1d_test_case_mod

  implicit none

  integer, parameter :: ns = 1             ! Stencil width

  real, allocatable, dimension(:,:) :: rho ! Tracer density being advected at cell centers
  real, allocatable, dimension(:  ) :: flx ! Flux at cell interfaces

  real :: a = 0.5d0                        ! Alpha parameter when it is 0.5, flux is Lax-Wendroff
  real c                                   ! CFL number, u * dt / dx
  character(10) :: flux_limiter_type = 'van_leer'
  integer i
  character(30) scheme

  namelist /params/ nx, nt, dt, u, flux_limiter_type

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,2))
  allocate(flx(1:nx+1))

  scheme = 'lax_wendroff_tvd.' // trim(flux_limiter_type)

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, time_step, ns, nx, x, rho(:,old))

  ! Run integration.
  c = u * dt / dx
  write(*, *) time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call flux(rho(:,old), flx)
    do i = 1, nx
      rho(i,new) = rho(i,old) - dt / dx * (flx(i+1) - flx(i))
    end do
    call apply_bc(ns, nx, rho(:,new))
    call advance_time()
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    write(*, *) time_step, sum(rho(1:nx,old))
  end do

  deallocate(flx, rho)

  call adv_1d_test_case_final()

contains

  subroutine sine_ic(x, rho)

    real, intent(in) :: x(nx)
    real, intent(out) :: rho(1-ns:nx+ns)

    real, parameter :: pi = atan(1.0d0) * 4.0d0
    integer i

    do i = 1, nx
      rho(i) = sin(x(i) * 2 * pi)
    end do
    call apply_bc(ns, nx, rho)

  end subroutine sine_ic

  subroutine flux(rho, f)

    real, intent(in) :: rho(1-ns:nx+ns)
    real, intent(inout) :: f(1:nx+1)

    integer i

    do i = 1, nx
      f(i+1) = u * rho(i) + 0.5d0 * u * (1 - c) * (rho(i+1) - rho(i)) * flux_limiter(rho(i-1:i+1))
    end do
    call apply_bc(1, nx, f)

  end subroutine flux

  real function flux_limiter(rho) result(res)

    real, intent(in) :: rho(-1:1)

    real, parameter :: eps = 1.0e-6
    real r

    r = (rho(0) - rho(-1) + eps) / (rho(1) - rho(0) + eps)

    select case (flux_limiter_type)
    case ('none')
      res = 1
    case ('minmod')
      res = max(0.0, min(1.0, r))
    case ('superbee')
      res = max(0.0, min(1.0, 2 * r), min(2.0, r))
    case ('van_leer')
      res = (r + abs(r)) / (1 + abs(r))
    case ('mc')
      res = max(0.0, min(2 * r, (1 + r) / 2, 2.0))
    case default
      write(*, *) '[Error]: Invalid flux_limiter_type ' // trim(flux_limiter_type) // '!'
    end select

  end function flux_limiter

end program tvd_adv_1d_case
